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INTRODUCTION  AND  OBJECTIVES 


Practical  applications  of  computational  fluid  dynamics  within  the  foreseeable  future 
must  necessarily  utilize  some  form  of  turbulence  modeling.  Computations  based  on  the 
Reynolds-averaged  equations  of  motion  require  all  turbulence  transport  of  momentum  and 
energy  to  be  modeled;  while  large-eddy  simulations  require  the  subgrid-scale  turbulence 
to  be  modeled.  In  both  cases,  the  uncertainties  in  modeling  turbulence  within  the  viscous 
sublayer  constitute  a  major  weak  link  in  the  overall  numerical  computation. 

During  the  past  two  decades  a  great  deal  of  new  experimental  information  has  been 
assembled  on  the  physics  of  organized  eddy  structures  in  turbulent  flow,  especially  within 
the  viscous  sublayer  (see,  for  example,  the  recent  review  of  Cantwell,  1981).  Yet  it  has  not 
been  possible  thus  far  to  incorporate  this  body  of  physical  information  within  the  frame¬ 
work  of  Reynolds-averaged  turbulence  modeling.  The  reason  is  fundamental,  reflecting  a 
well-known  limitation  in  the  Reynolds-average  approach  which  begins  by  time  averaging 
the  dynamic  equations  of  motion.  In  this  initial  mathematical  step  important  physical 
aspects  of  organized  eddy  structures,  such  as  phase  relationships  and  coherent  structure 
dynamics,  are  obliterated  irreversibly.  Consequently,  some  totally  different  approach  is 
required  if  the  observed  physics  of  coherent  eddy  motions  are  to  be  incorporated  within 
the  framework  of  a  turbulence  model. 

Quite  recently  a  new  approach  has  been  developed  for  modeling  viscous  sublayer  tur¬ 
bulence  in  incompressible  flow  without  heat  transfer  (Chapman  and  Kuhn  (1981).  Their 
method  models  directly  the  essential  organized  eddy  structures  observed  in  experiments. 
The  principal  steps  in  this  “coherent-structure  modelling”  are:  first,  to  model  velocity 
boundary  conditions  at  the  outer  edge  of  the  viscous  sublayer,  then  to  compute  time- 
dependent  dynamics,  and  finally  to  time  average  computed  results.  Thus,  time  averaging 
is  the  last  operation  performed  on  computed  dynamics,  rather  than  the  first  operation 
performed  on  the  dynamic  equations.  This  initial  effort,  although  far  from  fully  devel¬ 
oped,  already  has  been  surprisingly  successful  in  modeling  some  important  characteristics 


of  viscous  sublayer  turbulence. 


The  over  all  objective  of  the  present  research  is  to  develop,  using  the  Navier-Stokes 
equations,  a  computational  model  of  viscous  sublayer  turbulence  applicable  to  flow  with 
heat  transfer  and  to  compressible  flow.  Specific  tasks  within  this  objective  are  to  utilize 
the  model  to  compute  the  distribution  of  turbulent  Prandtl  number  within  the  viscous 
sublayer  (a)  for  fluids  of  various  molecular  Prandtl  number,  and  (b)  for  flows  with  adverse, 
zero,  and  favorable  streamwise  pressure  gradients. 

Navier-Stokes  computations  of  the  turbulent  Prandtl  number  Prt  are  significant  be¬ 
cause  experimental  techniques  have  been  unable  to  provide  reliable  measurements  within 
the  viscous  sublayer.  This  parameter,  of  course,  affects  turbulent  heat  transfer.  Since  by 
definition, 


Pr,= 


0»>)(g) 
W  (f ) 


an  experimental  determination  would  require  accurate  measurements  near  a  wall  of  the 
Reynolds  stress  ut>,  the  correlation  Bv  between  temperature  6  and  normal  velocity,  the 
mean  velocity  gradient  dU /dy,  and  the  mean  temperature  gradient  dB/dy.  To  date,  it  has 
not  been  possible  to  make  reliable  measurements  of  these  four  quantities  in  the  viscous 
sublayer.  The  band  of  uncertainty  from  one  set  of  data  for  air  (Simpson  et  al.  (1970)), 
is  illustrated  in  fig.  1.  Other  sets  of  experimental  data,  e.g.  Pulachier  (1972),  indicate 
that  the  overall  uncertainty  band  is  even  broader  than  this  one  set  of  measurements  would 
indicate.  Experiments  are  unable  to  determine  whether  the  values  of  Pr(  near  a  wall  are 
large,  small,  or  intermediate.  Experiments  also  have  not  been  able  to  define  how  Prt  varies 
with  molecular  Pr  or  with  pressure  gradient  within  the  viscous  sublayer.  We  believe  that 
Navier-Stokes  computations  can  shed  much  light  on  these  uncertainties. 


MOTIVATIONS  FOR  COHERENT  STRUCTURE  MODELING  OF 
VISCOUS  SUBLAYER  TURBULENCE 


Coherent-structure  modeling  departs  markedly  from  Reynolds-averaging  modeling, 

N 

and  thereby  offers  the  potential  of  some  entirely  new  advances  in  turbulence  computation. 
Both  the  motivations  and  the  payoff  for  this  type  of  research  are  quite  different  from  those 
for  more  conventional  turbulence  modeling.  In  addition  to  the  specific  motivations  outlined 
above  for  the  objectives  of  the  present  research,  there  are  other  important  motivations  that 
will  affect  different  areas  of  turbulence  computation.  It  is  appropriate,  therefore,  to  outline 
some  of  these. 

Three  other  motivations  for  developing  a  realistic  coherent-structure  model  of  viscous 
sublayer  turbulence  are: 

1  -To  provide  a  basis  for  strengthening  a  major  weak  link  in  present  Reynolds-average 
closure  schemes.  Conventional  turbulence  closure  methods  (e.g.  k— e,  other  2-equation 
methods,  and  Reynolds-Stress  transport  methods)  all  employ  differential  transport 
equations  that  are  “modeled”  forms  of  exact  transport  equations.  The  exact  equation 
for  the  transport  of  dissipation,  for  example,  contains  very  complex  turbulence  terms 
involving  pressure-velocity  correlations,  triple  correlations  of  velocity  gradients,  cor¬ 
relations  of  second  derivatives  of  velocity,  and  correlations  between  pressure  gradient 
and  velocity  gradient.  Such  correlations  are  not  measurable  with  present  experimental 
technology.  Because  modeled  equations  for  free  turbulence  yield  demonstrably  incor¬ 
rect  results  near  a  wall,  various  ad  hoc  functions  (up  to  5  in  number)  are  added  in  an 
effort  to  mend  this  shortcoming.  Without  any  guide  or  test  basis  from  experiment,  the 
inevitable  consequence  is  that  different  models  with  different  ad  hoc  functions  yield 
different  results  for  technically  important  quantities  such  as  skin  friction  in  flows  with 
pressure  gradient  (Patel  et.  al  1981).  If,  however,  a  realistic  Navicr-  Stokes  computa¬ 
tional  model  were  developed  for  the  time  dependent  dynamics  of  viscous  sublayer  flow, 
then  all  the  terms  in  the  exact  transport  equations  could  be  computed.  When  compared 
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with  corresponding  modeled  terms,  a  rational  formulation  would  be  obtained  of  the  wall 
damping  functions  for  integrating  turbulence  models  across  the  viscous  sublayer.  This 
would  strengthen  a  weak  link  in  conventional  turbulence  models.  Reynolds-average 
modeling,  of  course,  will  be  the  mainstay  of  practical  turbulent  flow  computation  for 
years  to  come. 

2  -To  provide  a  simple  test  flow  against  which  subgrid  scale  models  of  turbulence  in  large 
eddy  simulations  can  be  tested.  The  simple  flow  of  homogeneous  shear,  for  example, 
computed  in  detail  from  the  time-dependent  Navier-Stokes  equations  (Rogallo  1977) 
has  been  used  effectively  as  a  test  flow  for  assessing  subgrid  scale  turbulence  models 
for  application  in  the  outer  region  of  a  turbulent  boundary  layer  or  in  free  shear  layers 
(e.g.,  Clark  et.  al.  1979,  McMillan  and  Ferziger  1979,  Shirani  1981).  A  realistic 
coherent-structure  model  of  viscous  sublayer  turbulence,  likewise  computed  from  the 
time-dependent  Navier-Stokes  equations,  could  similarly  be  used  to  test  subgrid-scale 
models  for  large  eddy  simulations  of  turbulent  flow  in  the  region  adjacent  to  a  wall. 

3  -  To  provide  a  guide  for  modeling  the  lower  boundary  conditions  for  the  outer  turbulent 
region  in  large  eddy  simulations  of  boundary  layer  flow.  If  a  realistic  boundary  con¬ 
dition  of  this  type  could  be  developed,  the  entire  viscous  sublayer  could  be  modeled 
rather  than  directly  computed.  The  required  computer  power  for  large  eddy  simula¬ 
tions  at  high  Reynolds  numbers  would  thereby  be  reduced  by  a  very  large  factor  (about 
103;  Chapman,  1980). 

In  summary  of  the  motivations  for  developing  coherent-structure  turbulence  models, 
it  is  clear  that  the  potential  payoffs  are  significant  and  of  much  broader  scope  than  the 
specific  objectives  of  the  present  research. 


FORMULATION  OF  COMPUTATIONAL  MODEL 

DIFFERENTIAL  EQUATIONS 

The  numerical  scheme  selected  to  solve  the  Navier-Stokes  equations  is  an  adaptation 
of  the  Pulliam  and  Steger  (1980)  code  for  compressible  flow.  Velocities  in  the  physical  as, 
y,  z  domain  are  u,  v,  W.  In  the  computational  £,  rf,  f,  domain  the  dependent  variables 
are  p,  pu ,  pv,  pw,  and  p  (e,  +  (u2  +  v2  +  tu2)  / 2],  where  e,-  is  the  internal  energy  per  unit 
mass.  Pulliam  and  Steger  use  z  as  the  coordinate  normal  to  the  surface,  and  y  as  the 
spanwise  coordinate.  In  order  to  minimize  code  changes,  we  use  this  coordinate  system  in 
the  numerical  computations  and  then,  at  the  end,  change  z  to  y,  w  to  v,  y  to  z,  and  v  to  tu 
in  order  to  present  results  in  the  more  familiar  boundary  layer  coordinate  designations.  We 
employ  simple  cartesian  coordinates  in  the  streamwise  x-  and  spanwise  z  =  y-directions, 
and  a  stretched  mesh  in  the  normal  y  =  5  direction  in  order  to  concentrate  points  both  near 
the  wall  and  the  outer  edge.  The  metrics  for  the  transformation  are  £,  =  1,  —  1/Ay+, 

=  1/Az+  where  Ay+  and  A2+  are  the  mesh  spacings  in  wall  units  in  the  spanwise  and 
normal  directions  respectively. 

The  same  basic  physical  approximation  is  made  as  in  Chapman  and  Kuhn  (1981), 
namely,  that  coherent  sublayer  eddy  structures  are  highly  elongate  streamwise.  In  our 
case  we  assume-in  accordance  with  experimental  observations  (e.g.  Iritani  et.  al  (1981)) 
—  that  the  temperature  eddies  are  also  highly  elongate  streamwise.  In  the  notation  of 
Pulliam  and  Steger,  the  corresponding  divergence  form  of  the  Navier-Stokes  equations  is 


d,q  =  dn(Fv  -  F)  +  df{Gv  -G)  +  Fp 


(1) 


A  A  A  A 

where  the  vectors  F,  F„,  G ,  Gu  and  q  are  identical  to  the  conventional  ones  used  by  Pulliam 
and  Steger  for  the  case  of  orthogonal  coordinates  corresponding  to  the  metrics  of  equation 


(1).  However,  the  vector 


0 

_?£i 

dz 


Fp  =  J~l 


_ Bp, 

dz 


+  £(*«)]  J 


(2) 


differs  from  theirs  due  to  our  assumption  of  highly  elongate  eddies  in  the  streamwise 
direction.  Thus,  streamwise  derivatives  other  than  the  mean  streamwise  pressure  gradient 
dpjdx  are  neglected  relative  to  derivatives  in  the  spanwise  and  normal  directions. 

This  basic  approximation,  which  involves  computation  of  3  velocity  components  in  2 
space  directions,  has  been  variously  termed  “2.5D  flow”,  or  “slender  turbulence  theory”, 
or  “slender  eddy  theory.” 

Since  Pulliam  and  Steger  used  the  thin-layer  form  of  the  Navier-Stokes  equations,  we 
write,  in  order  to  retain  a  tridiagonal  structure, 


d(Fv  -  F)  _  5F 
dq  ~  dq 


"b 


d{G„  -  G)  =  dG 
dq  ~  dq 


—  Ba  +  Bx 


(3) 


where  the  matrices  A0  and  B0  contain  no  cross  derivative  terms,  while  Ax  and  Bx  contain 
all  such  terms.  The  detailed  structure  of  these  matrices  is  given  in  Appendix  I.  Our  basic 
factored  algorithm  becomes 


(I  -  !«„<)(/  -  +  a,cr  I  +  2-hd„A"A," 


+  lh6<B;Ag"  +  hk  +  lAq" 


(4) 


This  differs  from  the  Pulliam-Steger  algorithm  by  the  addition  of  the  last  four  terms  on 
the  right  hand  side.  Of  these,  the  last  one  corresponds  to  a  three-point  backward  Euler 

difference  scheme  for  time  derivatives  (added  for  accuracy  and  stability  reasons  explained 

\ 

later);  the  next  to  last  term  corresponds  to  body  pressure  gradients  imposed  by  the  large 
scale  eddies,  while  the  third  and  fourth  from  last  represent  the  additional  terms  required 
for  the  full  Navier-Stokes  equations  in  place  of  the  thin-layer  equations  used  by  Pulliam 
and  Steger. 


BOUNDARY  CONDITIONS 

The  primary  element  in  formulating  our  computational  model  involves  the  construction 
of  time-  and  space-dependent  boundary  conditions  that  are  as  realistic  as  possible  for  the 
fluctuating  temperature  6  and  the  fluctuating  velocity  components  u,  v,  w  at  the  outer 
edge  of  the  viscous  sublayer.  Various  experimental  observations  are  used  as  a  guide  to  this 
construction.  Computations  using  the  time-dependent  Navier-Stokes  equations  are  then 
made  over  a  length  of  time  long  enough  for  the  various  statistical  quantities,  e.g.,  mean 
values,  correlation  coefficients,  rms  fluctuating  intensities,  etc.,  to  reach  a  periodic  state 
within  the  viscous  sublayer. 

Our  general  approach  is  an  extension  of  that  employed  by  Chapman  and  Kuhn  (1981) 
for  incompressible  flow  without  heat  transfer.  They  considered  two  coherent  components 
of  turbulent  motion  at  the  outer  edge  of  the  viscous  sublayer  (lower  edge  of  the  logarithmic 
region).  One  component  represents  the  small  scale  eddies  (SSE)  responsible  for  the  prin¬ 
cipal  production  of  turbulence  and  Reynolds  stress.  A  second  component  represents  the 
organized  large  scale  eddies  (LSE)  which  interact  with  the  small  scale  eddies.  They  used 
a  simple  initial  construction  for  purposes  of  illustrating  the  potential  and  main  character¬ 
istics  of  this  type  of  turbulence  modeling.  With  all  quantities  expressed  in  conventional 
dimensionless  wall  variables,  their  velocity  boundary  conditions  for  the  outer  edge  of  the 
viscous  sublayer  were 


8 


Component  1-SSE 


Component  2-LSE 


ue+  =  2at  sin  NtT sin  f 


+  v/2(o2  -  af)  sin(Nu2T  +  <f>u2) 


ve+  =  —2(3  sin  NtT  sin  f 


(5) 


iue+  =  2(3  cos  NiT  cos  f 


+  T  +  t.,) 


where  a,  /?,  7  are  the  rms  intensities  of  fluctuation  of  ut+,  ve,  twe,  respectively,  JVj  is  the 
mean  frequency  of  SSE  burst  events  (ejection/sweep  events),  Nu2  is  the  mean  frequency 
of  LSE,  f  =  2 txz/\  is  a  dimensionless  spanwise  distance  variable,  T  is  the  dimensionless 
time,  A  is  the  mean  spacing  between  high-speed  or  low-speed  streaks,  at/a  is  the  uv 
correlation  coefficient  (0.45),  and  <^u2  and  <f>w2  are  phase  angles.  The  flow  model  is  periodic 
both  in  time  and  span.  Hence  the  boundary  conditions  for  the  spanwise  side  walls  of  the 
computational  domain  were 


u«+(y,<M)  =  ue+(y,A,<) 


ve+(y>  0,0  =  ve+(y,A  ,t) 


(6) 


u>e+(y,o,t)  =  u>e+(y,A  ,t) 


and  at  the  surface 


u*+(0, z, t)  =  ve+(0,z,t)  -  twe+(0 ,z,t)  =  0 


This  construction  of  velocity  boundary  conditions  was  guided  by  experimental  obser¬ 
vations  of  organized  sublayer  structure.  For  example,  the  sin  f  factor  for  ut+  corresponds 
to  the  observation  of  high-  and  low-speed  streaks  spaced  spanwise  a  mean  distance  of 


A  apart,  while  the  factors  sinf,  sinf,  cosf  in  the  ue+,  t/e+,  toe+  equations,  respectively, 
correspond  to  a  simple  modeling  of  the  observations  of  contra-rotating  vortical  motion 
near  a  wall.  The  180  deg  phase  difference  between  ue+  and  ve+  corresponds  to  observations 
from  conditional  samples  in  the  sublayer  that  u  and  v  are  180  deg  out  of  phase  during  the 
Reynolds-stress  intensive  ejection-sweep  event.  The  90  deg  phase  difference  in  time  be¬ 
tween  vt  and  tue  corresponds  to  the  observation  that  the  derivative  ( dv2/dy+)t  at  the  outer 
edge  of  the  viscous  sublayer  is  zero.  A  comparison  of  a  number  of  turbulence  characteristics 
computed  from  this  simple  model  with  those  measured  in  experiments  showed  surprisingly 
good  agreement.  This  model  was  developed  for  flows  with  a  mean  streamwise  pressure 
gradient,  although  numerical  computations  were  made  by  Chapman  and  Kuhn  only  for 
the  relatively  small  pressure  gradients  that  exist  in  incompressible  pipe  and  channel  flow. 

The  above  two-component  velocity  boundary  conditions  at  the  outer  edge  of  the  vis¬ 
cous  sublayer  are  fully  compatible  with  the  concept  of  “active”  and  “inactive”  components 
of  turbulent  motion  in  the  log  region  as  characterized  by  Townsend  (1961)  and  Bradshaw 
(1967).  Component  1  is  active,  involving  small  scale  eddies,  producing  the  Reynolds  stress, 
being  rotational,  and  being  dependent  on  wall  variables.  Component  2  is  inactive,  involv¬ 
ing  large-scale  eddies,  producing  energy  but  no  Reynolds  stress,  being  irrotational,  and 
being  dependent  on  outer  variables.  There  is  much  experimental  information  that  turbu¬ 
lent  flow  in  the  log  region,  and  hence  on  the  outer  edge  of  the  viscous  sublayer,  comprises 
these  two  distinct  types  of  motion. 

The  computational  model  of  Chapman  and  Kuhn  is  not  regarded  as  sufficiently  real¬ 
istic  for  purposes  such  as  outlined  in  the  section  above  on  motivations.  Their  results  do 
illustrate,  however,  that  the  general  approach  employed  has  the  potential  of  eventually 
leading  to  a  sufficiently  realistic  model.  Towards  this  end,  it  is  believed  that  additional 
aspects  of  turbulence  physics  will  have  to  be  introduced  into  the  model.  For  example, 
production  of  Reynolds  stress  should  be  intermittent  in  time,  rather  than  sinusoidal-like 
as  in  boundary  conditions  (5).  Also,  conditional  samples  for  v[t)  in  ejection/sweep  events 
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indicate  that  the  v 2  turbulent  energy  contained  in  these  SSE  events  should  be  only  a  frac¬ 
tion  of  the  total  v 2  energy,  rather  than  the  entire  amount  as  in  equation  (5).  Further, 
visualizations  of  viscous  sublayer  flow  exhibit  a  mixture  of  order  and  disorder,  rather  than 
the  coherent  order  only  of  equation  (5).  For  the  heat  transfer  case  we  must  construct 
appropriate  correspondingly  boundary  condition  for  the  fluctuating  temperature  9e  at  the 
outer  edge  of  the  viscous  sublayer. 

One  important  guide  to  the  construction  of  boundary  conditions  is  provided  by  mea¬ 
surements  of  spectral  density  /.  Fortunately,  spectra  for  0,  u,  v ,  and  w,  have  been  obtained 
by  Fulachier  (1972)  at  y-positions  near  the  outer  edge  of  the  viscous  sublayer.  Data  for 
ye+  «  40  (interpolated  between  his  measurements  at  y+  =  31  and  77)  are  shown  in  figure 
2.  Since 

r  fdk  =  r  Jb/d(ln  Jfc)  =  1 

Jo  Jo 

the  areas  under  curves  of  u2kfu,  v2kfV)  and  tv2kf„  versus  log  A:  are  proportional  to  the 
relative  amounts  of  kinetic  energy  in  these  velocity  components.  Along  the  wave  number 
scale  k,  pips  are  shown  corresponding  to  LSE  eddy  scales  of  tt/<5  (where  6  is  the  boundary 
layer  thickness)  and  to  SSE  eddy  scales  of  tt/A.  It  is  seen  that,  whereas  both  u  and  9 
exhibit  eddy  scales  ranging  from  large  [0(5)]  to  small  [0(A)],  v  exhibit  scales  ranging  only 
from  medium  [O  (several  A)]  to  small.  Our  boundary  conditions  will  be  constructed  to  be 
as  consistent  as  possible  with  experimental  data  of  this  type. 

A  relatively  simple  method  (developed  under  a  separate  contract)  is  used  to  simulate 
the  intermittent  character  of  Reynolds-stress  production.  Since  uv  is  produced  only  by 
the  active  component  of  small  scale  eddies  (SSE),  only  the  time  functions  sinATjT  and 
cos  NtT  in  equation  (1)  need  be  modified.  For  this  purpose  we  use  a  simple  Fourier  series 
approximation  //(<£)  to  a  rectangular  pulse  function,  where  (f>  =  NiT. 
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For  M  terms  in  the  approximating  Fourier  series,  we  have 


values  of  M  in  the  range  of  3  to  5  are  anticipated.  Initial  computations  will  be  made  with 
Xi  =  X2  =  0.25jt,  since  hot-wire  measurements  have  indicated  that  the  principal  ui)  is 
produced  in  about  25%  of  the  time,  with  relative  uv  quiescence  the  remaining  75%  of  time. 
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In  conventional  wall  variables,  the  outer  edge  boundary  conditions  constructed  to  date 
include  terms  representing  small  scale  eddies  (SSE)  of  length  scale  A,  large  scale  eddies 
(LSE)  of  length  scale  >  10A,  and  medium  scale  eddies  (MSE)  of  intermediate  length  scale 
3A.  Computations  are  made  within  a  domain  covering  3  spanwise  cells  each  A  in  width. 
The  LSE  are  treated  as  being  of  sufficiently  large  scale  that  their  time-dependent  velocity 
component  does  not  vary  across  the  spanwise  extent  of  the  computational  domain.  The 
equation  of  continuity  precludes  having  a  LSE  component  in  the  boundary  condition  for 
vt. 

SSE  LSE  MSE 

Scale  A  Scale  >10A  Scale  3A 

ue  =  \/2a,,Fu(^)sin^  +  \j2{o?  -  a?)  sin (Nu2T  +  <f>u2) 

ve  =  \Z2/?i Fv(<f>)  sin  $■  +  2p2  sin  (3Nu2T)  sin  | 

o 

we  =  \[2 PiFw((f))  cos  f  -  Pi)  sin 

(8) 

9t  =  \fl alFe{4>)  sin  ?  +  \[l{a 2  “  °i)  sin  {Nu2T  +  <f> 02 )  (9) 

where,  in  addition  to  the  quantities  already  defined,  a  =  \f[02)t ,  <j>o2  is  a  phase  angle  in 
the  large-eddy  component  of  temperature,  and 


m) = -km = F.M  =  "M 


r»  / 


HU  +  #„,) 


(10) 


are  the  intermittent  time  functions  normalized  such  that  F2  =  1.  The  boundary  conditions 
for  the  spanwise  side  walls  of  the  computational  domain  are 

u«+(y,0  ,t)  =  u«+(y,3A,i) 
w«+(y.o,0  =  ue+(y,3A,<) 
t"«+(y>M)  =  u>e+(y,3A  ,t) 

Thus  the  computation  allows  for  flow  interaction  between  the  center  cell  and  the  two 
outside  cells  unconstrained  by  specific  boundary  conditions. 

As  in  the  model  of  Chapman  and  Kuhn:  a,  /?,  7,  a  are  determined  from  experimental 
data  on  the  intensities  of  turbulence  fluctuations  in  wall  units  at  the  outer  edge  of  the 
viscous  sublayer  (a  «  2,/3  «  1,7  «  1.3, a  «  1.3);  <£tt2,  $02,  <Pu2,  are  determined  by 
computer  trial  to  conform  as  well  as  possible  to  the  mean-velocity  law  of  the  wall,  mean- 
temperature  law  of  the  wall,  and  u-velocity  skewness,  respectively;  </vi  is  determined 
analytically  by  the  requirement  that  ( dv2jdy)t  is  zero;  07  is  determined  by  the  correlation 
coefficient  0.45  =  —  (Ru„)e  =  aipl/aP\  aj  is  determined  by  the  correlation  coefficient 
0.8  «  -{Rll9)e  =  “o"1  +  oon/O*2  ~  «i)(a2  -  o.\)  cos  (^92  -  ^2);  and  PJP  is  determined 
from  conditional  samples  of  the  ejection/sweep  event.  This  latter  determination  is  made 
by  assuming  that  the  peak  to  peak  amplitude  ratio  <  v  >ppa  /  <  u  >ppa,  from  the 
conditional  samples  is  equal  to  the  ratio  Pi/alt  and  then  using  the  above  equation  for 
(*■■).  to  compute  both  Pi/P  and  ak/a.  An  independent  determination  of  at/a  has  been 
made  from  direct  measurements  of  the  fraction  of  total  u2  energy  that  exists  during  bursts, 
and  hence  also  of  PJP  from  the  equation  for  (Rvv)e.  Still  a  different  determination  has 
been  made  from  <  w  >pp0  /  <  u  >pp(I=  Pi/oti  coupled  with  the  equation  for  (f2„t,)e. 


Results  of  these  determinations  (made  on  a  separate  contract)  are: 


Data  Source 

Method 

fii/fi 

Chen  Ac  Balckwelder  (1978) 

<v>  /  <u>  plus  R„v  =  -.45 

.53 

Nakagawa  Ac  Nezu  (1981) 

<  v  >  /  <  u  >  plus  =  -.45 

.72 

Blackwelder  &  Kaplan  (1979) 

<  v  >  (  <  u  >  plus  R,,,,  =  -.45 

.49 

Kim  (1982) 

<  u  >  /  <  v  >  from  LES  computations, 
plus  RuV  =  -.45 

.60 

Kim,  Kline  Ac  Reynolds 
(1971) 

Fraction  of  total  u2  energy  during 
bursts  (.68),  plus  R„v  =  —.45 

.59 

Blackwelder  Ac  Kaplan  (1979) 

<  w  >  /  <  u  >  plus  Ru„  =  —.45 

.64 

On  average,  the  value  of  appears  to  be  near  0.6. 

It  is  to  be  noted  that  the  structure  of  boundary  conditions  (8)  comprises  in 
equal  proportions  both  ejection/sweep  events  (where  a  sweep  follows  an  ejection),  and 
sweep/ejection  events  (where  an  ejection  follows  a  sweep).  The  recent  experiments  of 
Johansson  and  Alfredsson  (1982)  clearly  reveal  that  both  types  of  event  occur.  Their  ex¬ 
perimental  VITA  data  for  the  smallest  threshold  values  and  the  longest  integration  times 
employed,  indicate  these  two  event  types  to  be  about  equally  numerous.  Both  types  have 
also  been  observed  visually  by  Offen  and  Kline  (1974)  who  use  the  term  “cleansing  sweep” 
in  an  ejection/sweep  event  to  distinguish  from  sweeps  in  a  sweep/ejection  event. 


STATUS  OF  RESEACH 


I.  Recapitulation  of  Research  During  1982 

During  the  first  year  of  research,  effort  was  directed  mainly  on  surveying  various  ex¬ 
isting  theories  for  Prt  in  the  viscous  sublayer,  on  programming  a  Navier-Stokes  code,  and 
on  testing  and  debugging  this  code. 

Survey  of  various  theories  for  Prt 

Without  a  reliable  guide  from  the  experiment  as  to  how  Prt  should  be  modeled  near 
a  wall,  the  various  theoretical  models  advanced  to  date  have  differed  greatly.  This  is 
illustrated  in  figures  3,  4,  and  5  showing  a  compilation  of  various  theories  for  fluids  of 
molecular  Prandtl  number  0.72  (air),  0  (water),  and  1000  (oil),  respectively.  These  various 
theories  are  for  zero  pressure  gradient  and  clearly  demonstrate  that  a  very  wide  range  of 
uncertainty  exists  for  Prt  in  the  viscous  sublayer.  Present  theories  about  the  dependence 
of  Pr,  on  pressure  gradient  are  fewer  in  number  but  also  differ  greatly. 

Programming  of  Navier-Stokes  Code 

The  implicit,  compressible-flow,  3-D  code  of  Pulliam  and  Steger  (1980)  has  been 
used  as  the  basic  Navier-Stokes  code  to  modify  and  adapt  for  viscous  sublayer  modeling. 
This  code  has  been  modified  in  five  ways  as  follows:  (1)  the  grid  has  been  stretched  to 
cluster  points  near  the  wall  and  the  outer  edge;  (2)  the  flow  variables  have  been  made 
nondimensional  in  the  conventional  wall  variables  used  for  turbulent  flows;  (3)  various 
necessary  viscous  terms  have  been  added  to  bring  the  Pulliam-Steger  “thin-layer”  code 
into  the  full  Navier-Stokes  form;  (4)  the  boundary  conditions  have  been  restructured  to  be 
appropriate  for  computations  of  viscous  sublayer  turbulence;  (5)  the  algorithm  for  time 
derivatives  has  been  revised  for  increased  accuracy. 

Testing  and  Debugging  of  Code 

In  order  to  test  the  algorithm  accuracy  for  the  time  derivative  terms,  an  exact  an¬ 
alytical  solution  was  derived  for  incompressible  oscillating  shear  flow  with  heat  transfer. 
This  solution  is  an  extension  of  the  Chapman-Kuhn  solution  for  flow  without  heat  trans- 
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fer.  Mathematical  details  of  the  analytical  solution  are  presented  in  Appendix  II.  It  was 
found  that  neither  of  the  option  algorithms  in  the  original  Pulliam-Steger  code  (i.e.,  Euler 
implicit  first  order,  or  trapezoidal  second  order)  were  sufficiently  accurate  for  turbulence 
computations.  With  helpful  guidance  from  Dr.  Pulliam,  a  three-point  backward  Euler 
algorithm  was  coded  and  found  to  be  satisfactory.  The  very  close  agreement  obtained 
with  this  second-order  accurate  algorithm,  between  numerical  computations  and  the  ex¬ 
act  analytical  solution,  is  illustrated  in  figures  6  and  7  for  the  adiabatic-wall  case  and 
the  heat-transfer  case,  respectively.  Since  oscillating  shear  flow  with  heat  transfer  involves 
viscous  and  heat-conduction  terms  as  well  as  time-derivative  terms,  this  test  flow  provides 
a  check  on  the  accuracy  of  the  numeries  for  all  three  types  of  terms.  It  does  not  provide, 
however,  a  check  on  the  accuracy  of  the  important  non-linear  convection  terms  which  are 
zero  in  oscillating  shear  flow. 

Our  initial  runs  of  the  viscous  sublayer  code  revealed  an  instability  problem  associated 
with  a  large  residual  in  the  continuity  equation,  and  a  “decoupling”  of  the  computed  data 
in  the  y-direction  normal  to  the  wall.  In  this  decoupling,  computed  values  of  v,p  and  p' 
oscillated  from  grid  point  to  grid  point.  Decoupling  is  usually  a  signal  of  a  problem  with 
boundary  conditions.  (This  instability  was  minimized  to  acceptable  levels  by  a  coding 
change  described  below  in  the  section  on  progress  during  1983.)  Examination  of  the 
continuity  residual  revealed  that  the  instability  arose  at  the  outer  edge  of  the  viscous 
sublayer,  since  there  the  residual  was  around  10-1,  whereas  in  the  inner  main  portion  of  the 
flow  it  was  the  order  of  10-3  to  10-c.  This  also  suggests  that  the  initial  boundary  conditions 
at  the  outer  edge  were  creating  a  problem.  Such  problems  can  arise  form  either  physical  or 
numerical  causes.  Some  trials  with  obviously  wrong  boundary  conditions  produced  very 
extreme  decoupling,  leading  to  a  blow  up  of  the  code. 

II.  Research  Progress  During  1983 

During  the  past  year,  research  has  been  focused  on  completing  the  programming,  on 
investigating  various  types  of  boundary  conditions,  and  on  conducting  some  preliminary 
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computer  runs  to  explore  the  effects  on  turbulent  Prandtl  number  of  variations  in  molecular 
Prandtl  number  and  in  mean  streamwise  pressure  gradient. 

Programming  Modifications 

Considerable  code  modification  was  undertaken  in  order  to  expand  the  data  output 
and  graphics  program,  to  increase  code  robustness,  and  to  transfer  the  code  to  two  new 
computers  installed  at  the  Ames  Research  Center  during  1983.  The  data  output  and  associ¬ 
ated  graphics  program  was  extended  to  handle  26  statistical  quantities:  the  mean  velocities 
u,  v,  tZ>,  temperature  0,  pressure  p  and  density  f ,  the  Reynolds  stresses  uv,  uw  and  velocity- 
temperature  correlations  u9 ,  v9 ,  the  intensities  of  fluctuations  u',  v\  w',  0',  p',  (uu)', 

the  correlation  coefficients  Ruv,  Rv0i  R^g,  the  skewness  Su,  Se ,  Suv  and  the  flatness  fac¬ 
tors  Fu,  Fg,  Fuv.  It  was  found  that  refinement  of  the  grid  from  25  xl7  to  39  x  25  in 
the  spanwise  plane  increased  considerably  the  codes’  robustness  in  handling  a  variety  of 
boundary  conditions  that  were  explored.  The  finer  grid  apparently  localizes  boundary- 
condition  problems  near  the  outer  edge  of  the  viscous  sublayer,  and  allows  for  a  greater 
portion  of  the  interior  flow  to  be  well  behaved.  Some  tests  were  made  of  adding  an  artificial 
dissipation  subroutine  to  the  code  in  order  to  assess  its  effects  on  the  stability  and  decou¬ 
pling  problems.  Both  stability  and  decoupling  were  improved,  but  the  numerical  results 
were  altered  considerably,  and  hence  the  artificial  dissipation  subroutine  was  dropped  from 
further  consideration. 

During  1983  NASA  Ames  replaced  the  CDC  7600  computer,  which  we  had  originally 
programmed  for,  first  by  the  CRAY--1S,  and  then  later  replaced  this  computer  by  the 
CRAY-XMP.  The  corresponding  changes  in  operating  system  required  some  modifications 
to  the  code.  To  date  all  three  of  these  computers  have  been  used  in  the  present  research 
project. 

Research  on  Boundary  Conditions 

In  a  code  for  3-D  compressible-viscous  flow,  five  boundary  conditions  are  required 
at  every  boundary.  Not  any  combination  of  five  conditions  will  work.  This  problem  is 
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particularly  acute  at  the  out  edge  of  our  computational  domain  where  non-zero  velocity 
components  are  specified  and  where  fluid  flows  in  and  out  across  this  boundary.  The  set 
of  boundary  conditions  which  is  physically  and  mathematically  correct  for  this  type  of  vis¬ 
cous  compressible  flow  is  not  known  a  priori.  Since  no  sound  theoretical  guide  is  available, 
considerable  computer  trial  is  necessary  using  various  types  of  boundary  conditions.  Fur¬ 
ther  complexity  arises  because  the  numerical  boundary  conditions  can  be  treated  either 
explicitly  or  implicitly.  To  date  explicit  boundary  condition  schemes  have  been  mainly 
used.  One  implicit  scheme  has  been  tried,  but  has  not  provided  an  improvement  in  either 
accuracy  or  stability. 

In  the  present  research  we  have  generalized  the  boundary  conditions  of  Chapman 
and  Kuhn  (1981)  who  worked  with  incompressible  flow.  Difficult  problems  of  boundary 
conditions  arise  when  it  is  attempted  to  compute  a  compressible  flow.  If  the  dissipative 
terms  were  to  dominate  at  the  outer  boundary,  then  a  stable  scheme  would  likely  occur 
with  the  specification  of  velocities  ue,  ve,  we,  temperature  0e,  and  either  pressure  pe  or 
density  pe.  If  the  non-linear  convective  terms  dominate,  however,  then  this  may  only  be 
weakly  stable  at  best.  Under  such  circumstances  a  characteristics  type  set  of  boundary 
conditions  might  have  to  be  implemented.  With  Navier-Stokes  codes  the  situation  is  not 
usually  near  either  of  these  two  extremes,  but  is  somewhere  in  between.  Specification  of 
wrong  boundary  conditions  can  result  in  serious  stability  problems.  In  the  present  research 
we  have  attempted  to  structure  boundary  conditions  that  use  as  much  as  possible  of  the 
previous  work  of  Chapman  and  Kuhn,  and  still  provide  adequate  stability. 

One  change  was  found  to  be  necessary  to  achieve  stability.  At  the  outer  edge,  the 
mass  flux  quantities  (pu)n  ( pw)e,  (pw)e  are  specified  instead  of  the  velocity  components. 
In  this  way  mass  is  conserved  within  the  computational  domain.  For  the  fourth  boundary 
condition  the  outer-edge  temperature  9t  has  been  modeled  in  a  structural  form  similar  to 
that  of  the  u-componeut  of  velocity.  For  the  fifth  boundary  condition  at  the  outer  edge  of 
the  viscous  sublayer,  there  are  a  number  of  possibilities  of  which  the  following  have  been 


investigated: 

(a)  pe  extrapolated  from  y-momentum  equation 

(b)  pe  extrapolated  from  continuity  equation 

(c)  pe  =  pe  (1  +  esm{NiT  +  ^i)  +  6  sin^T  +  ^2)1 

(d)  pe  =  pe  [1  +  esin (N{T  +  ^i)  +  6  sin {NtT  +  fa)]'1 

(e)  Pe  extrapolated  from  (|^)e  =  0 

It  is  noted  that  boundary  conditions  of  the  types  (a)  and  (b)  do  not  introduce  any  addi¬ 
tional  physical  assumptions  about  the  turbulence,  and  hence  are  the  most  desirable.  Thus 
far,  the  most  success  has  been  obtained  with  conditions  (a),  and  the  numerical  results 
presented  later  were  obtained  using  this  condition. 

Two  quite  different  types  of  boundary  condition  also  have  been  explored  wherein 
derivatives  of  velocity  are  modeled: 

(f)  (pv)e,  {ptv)e,  8t  specified  with  condition  (a)  and  (|“)e  =  (fpe(l  +  esm(Ari!T+0l)  sin  £  + 
Ssin(NtT  +  fa)] 

(g)  ( Pu)c ,  (Hr,  (p«0 r  specified  with  condition  (a)  and  (g)e  =  (g|)e[l  +  esin^T  + 
<f>i)  sin  £  +  6  sin(7ViT  +  <f>2 )] 

These  two  types  of  boundary  conditions  yielded  somewhat  more  realistic  values  for  the 
temperature  intensity  6'  and  the  fluctuating  pressure  intensity  p'.  Other  quantities  such 
as  mean  velocity  u  and  mean  temperature  9,  however,  were  less  realistic. 

One  further  generalization  in  the  structure  of  the  outer-edge  boundary  conditions  was 
coded  in  1983.  The  Fourier  series  equations  representing  intermittent  Reynolds  stress  pro¬ 
duction,  equations  (7)  and  (10),  have  been  included  in  the  computer  program.  Exploratory 
cases  with  this  type  of  boundary  condition  will  be  run  during  1984. 

Preliminary  Results  for  Various  Pressure  Gradients  and  Molecular  Prandtl  Number 
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as 

0«  =  2cq  sin(JVxT  +  <f>n)  sin  £  +  \f2(a2  -  a])  s\n(Nu2T  +  <f>i2) 

Here  o  is  the  intensity  of  temperature  fluctuation  at  the  outer  edge  (a  =  1.5  in  wall 
variables).  The  phase  angles  <j>glt  <f>g2,  and  the  coefficient  ax  for  these  preliminary  results 
were  selected  so  that  the  u9  correlation  coefficient  would  be  —0.8  at  the  outer  edge  and 
the  v9  correlation  coefficient  0.45.  These  preliminary  results  serve  to  illustrate  the  type  of 
results  to  be  produced  by  our  Navier-Stokes  code,  as  well  as  the  qualitative  dependence 
of  Pr,  on  the  parameters  such  as  pressure  gradient  and  molecular  Pr. 

In  figure  8  a  comparison  is  presented  of  our  preliminary  computations  with  experimen¬ 
tal  results  for  mean  temperature  0,  intensity  of  temperature  fluctuation  O',  and  the  two 
temperature-velocity  correlation  coefficients  Rug  and  Rv8.  These  data  are  for  air,  Pr  =  0.7, 
with  zero  pressure  gradient  and  an  isothermal  wall.  The  computed  mean  temperature  field 
and  Rug  agree  well  with  experiment,  but  the  computed  temperature  intensity  and  Rvg  are 
higher  than  experiment.  Computed  values  of  Rl.g  are  believed  to  be  high  because  all  of  the 
modeled  v  energy  in  the  boundary  conditions  of  equation  (5)  is  contained  in  the  small  scale 
eddy  component  which  is  correlated  with  the  small  scale  eddy  component  of  temperature. 
It  is  anticipated  that  use  of  a  more  realistic  boundary  condition  for  ve,  such  as  represented 
by  equation  (8),  wherein  only  part  of  the  v  energy  is  contained  in  the  correlated  small 
scale  component,  will  lower  the  computed  Rt,g  to  values  in  closer  agreement  with  experi¬ 
ment.  Reasons  for  the  relatively  high  values  computed  for  temperature  intensity  are  not 
understood  at  the  present  time.  Further  research  with  varied  boundary  conditions  should 
clarify  this. 

Some  computations  have  been  made  in  order  to  provide  a  preliminary  assessment  of  the 
effects  of  molecular  Prandtl  number  Pr  on  the  turbulent  Prandtl  number  Pr,  in  the  viscous 
sublayer.  Runs  were  conducted  for  Pr  =  0.7,  1.5,  3,  6.  Higher  values  of  Pr  encountered  a 
stability  problem  which  is  being  investigated.  The  computed  results  for  an  isothermal  wall 


21 


>k\v 


xiz 


1  l.'.J.,.|.W-,.V.'. 


are  presented  in  figure  9.  In  some  limited  respects  these  results  agree  qualitatively  with 
two  of  the  various  theories  advanced  to  date,  namely,  the  theories  of  Cebeci  (1970)  and 
Geshev  (1978).  The  summary  of  various  theories  in  figures  3,  4,  and  5  include  these  two. 
It  is  interesting  that  the  computed  effect  of  Pr  on  Prt  near  the  wall  is  opposite  to  that 
near  the  outer  edge  of  the  viscous  sublayer.  The  computed  curves  intersect  near  y+  —  5. 
Near  the  wall,  our  computed  Pr,  increases  considerably  with  Pr.  Geshev’s  theory  for  an 
isothermal  wall  indicated  a  similar  increase. 


Values  of  Prt  at  y+  =  2 
Present 


Pr 

Computation 

Geshev 

Cebeci 

0.7 

.98 

.85 

1.25 

6.0 

2.0 

2.2 

1.04 

Other  theories  among  those  compiled  in  figures  3-5,  however,  do  not  indicate  such  a 
strong  increase  of  Prt  with  increasing  Pr  near  a  wall;  and  Cebeci’s  theoiy  indicates  that 
Prt  near  the  wall  decreases  with  increasing  Pr.  It  is  noted  that  Geshev’s  treatment  of 
wall  temperature  boundary  condition  is  more  realistic  than  the  others,  and  affects  v9  and 
hence  Pr,  near  the  wall  in  an  important  way.  Near  the  outer  edge  of  the  viscous  sublayer 
the  effect  of  Pr  on  Pr,  is  relatively  small,  as  might  be  anticipated.  The  computed  trend 
of  Pr,  >  1  for  Pr  <  1  and  Pr,  <  1  for  Pr  >  1  is  consistent  with  experimental  data.  In 
this  case  Cebeci’s  theory  indicates  the  same  trend  but  Geshev’s  does  not. 

Values  of  Prt  at  y+  =  20 


Present 

Pr 

Computation 

Geshev 

Cebeci 

0.7 

1.12 

.95 

1.16 

6.0 

.85 

1.10 

1.00 

The  various  theories  do  not  exhibit  the  near  independence  of  Pr,  on  Pr  around  y+  «  5 
that  our  preliminary  computations  exhibit. 

Some  preliminary  computer  runs  were  also  made  to  explore  the  effects  of  mean  stream- 
wise  pressure  gradient  on  Pr,  for  air  (Pr  =  0.7)  flowing  over  an  isothermal  wall.  In  this 
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case  the  boundary  conditions  were  the  same  as  above,  except  that  the  bursting  frequency 
TVj  of  small  scale  eddies  was  changed  with  pressure  gradient  in  accordance  with  the  ex¬ 
perimental  data  of  Schraub  and  Kline  (1965).  The  values  selected  for  pressure  gradient 
in  wall  variables  were  dp+/dx+  =  —.02  (favorable),  0,  +.02  (adverse).  As  illustrated  in 
figure  10,  the  computed  results  indicate  relatively  little  effect  of  mean  streaniwise  pressure 
gradient  on  Prt.  In  this  case  also  the  qualitative  effect  of  dp/dx  on  Prt  near  the  wall  is 
opposite  to  that  near  the  outer  edge.  Near  the  outer  edge,  the  trend  of  increasing  Prt 
with  increasing  favorable  pressure  gradient  (decreasing  dp/dx  )  is  qualitatively  consistent 
with  the  trend  deduced  from  experimental  data  (Blackwell  (1972),  Gibson  (1981)).  This 
trend  is  the  opposite  from  that  indicated  by  the  recent  theory  of  Thomas  (1982).  Below 
is  a  comparison  of  Prt  values  near  the  wall. 

Values  of  Prt  at  y+  =  5 


dp-k 

Present 

Thomas 

Blackwell 

dx_+ 

Computation 

Theory 

Experiment 

-.02 

.95 

0 

1.12 

.85 

1.80 

+.02 

1.12 

1.00 

1.60 

Near  the  outer  edge  of  the  viscous  sublayer  the  corresponding  results  are: 

Values  of  Prt  at  y+  —  20 


dp+ 

Present 

Thomas 

Blackwell 

dx+ 

Computation 

Theory 

Experiment 

-.02 

1.20 

0 

1.12 

.90 

.85 

+.02 

1.00 

1.10 

.75 

Unlike  our  computed  results,  neither  the  results  of  Thomas  or  of  Blackwell  exhibit  the 
crossover  and  nearly  common  intersection  near  y+  «  12.5  of  the  curves  for  various  /dx+. 

It  is  emphasized  that,  since  our  present  quantitative  results  are  preliminary  and  will 
likely  change  somewhat  as  the  temperature  and  velocity  boundary  conditions  are  refined, 
the  precise  numerical  values  are  not  regarded  as  significant.  The  observed  qualitative 
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trends,  however,  such  as  the  strong  increase  near  the  wall  of  Prt  with  increasing  Pr,  and 
the  relatively  small  effect  of  dp/dx  on  Prt,  may  not  be  markedly  affected  by  such  changes 
in  these  boundary  conditions.  These  preliminary  results  serve  to  reveal  mainly  the  types 
of  information  produced  by  the  Navier-Stokes  code  being  developed,  and  the  qualitative 
trends  which  may  not  be  altered  as  boundary  conditions  on  temperature  and  velocity  are 
further  refined.  If  the  preliminary  result  of  only  small  effects  of  dpfdx  on  Prt  holds  up 
with  more  refined  boundary  conditions,  this  would  be  a  simplifying  and  helpful  result  for 
practical  computations  of  heat  transfer.  These  preliminary  results  serve  to  illustrate  the 
inconsistencies  and  deficiencies  of  the  various  theoretical  models  for  Prt  in  the  viscous 
sublayer. 


III.  Outline  of  Research  Anticipated  During  1984 

Details  of  the  research  effort  planned  during  1984  have  been  given  in  our  recent  pro¬ 
posal  dated  July  1983  for  continuation  of  the  present  research  project.  Very  briefly,  the 
main  effort  will  be  expended  in  the  areas  of:  (1)  refining  the  temperature  boundary  con¬ 
dition;  (2)  using  more  complex  velocity  boundary  conditions  such  as  those  representing 
intermittent  Reynolds  stress  productions;  and  those  involving  two  components  of  v  veloc¬ 
ity;  and  (3),  to  the  extent  that  time  permits,  investigating  ways  of  further  reducing  the 
decoupling  observed  in  v,  p,  and  p'  very  near  the  wall,  and  ways  of  further  improving  the 
stability  and  robustness  of  the  code. 
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APPENDIX  I 


DECOMPOSITION  MATRICES 


The  derivatives  dF/dq  and  dG/dq  are  each  decomposed  into  two  matrices,  one  (subscript  zero) 
that  does  not  contain  cross  derivative  terms,  and  another  (subscript  2)  that  contains  all  such  terms. 
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where  the  matrix  A0  is 
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and  the  matrix  Ax  is 
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APPENDIX  H 


An  exact  analytical  solution  has  been  derived  for  oscillating  shear  flow  with  heat  trans¬ 
fer  in  an  incompressible  fluid  with  constant  properties.  Since  the  momentum  equations 
are  independent  of  the  thermal  energy  equation  in  incompressible  flow,  the  velocity  field 
is  the  same  as  for  the  case  without  heat  transfer  derived  by  Chapman  and  Kuhn  (1981). 
With  y  as  the  coordinate  normal  to  the  infinite  flat  surface,  and  u  and  v  as  the  velocity 
components  in  the  streamwise  and  spanwise  directions,  respectively,  we  have 

u  —  A  [sin  nt  -  e“*ntf  sin (nt  -  +  cy 

w  =  B  [sin(mt  +  ( <f> )  -  e~kmV  sin(mf  -  fcmy)] 


where 


k 


II 


The  oscillating  frequencies  are  n  in  the  streamwise  direction  and  m  in  the  spanwise  direc¬ 
tion,  and  the  mean  velocity  profile  is  u  =  cy.  For  this  velocity  field,  the  solution  to  the 
thermal  energy  equation 


pcp 


d6 

dt 


is 


T(z,  t)  =  Biy 2  +  clV  +  c2  +  B;2k^  +  Bze~2k ™* 

+  A2e~knV  [sin(n<  -  kny)  -  cos(nt  -  &„y)]  -  Aie~2knV cos2(nt  -  kny) 
-  A4e~2krnV  cos2(mf  +  <f>  -  kmy)  -  A2e~9lV  sin(nf  -  pty) 

—  A2e~9>l  cos  {nt  -  p4y)  +  A3e~9iV  cos(2nf  -  p2y) 

+  i4<e“^8|,cos(2mf  +  2 <f>-  p3y) 
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Figure  5.  Some  theoretical  models  of  turbulent  Prandtl  number 
for  a  molecular  Prandtl  number  of  1000  (oil) . 
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;ure  1.  Comparison  of  numerical  computations  with 
exact  analytical  results  for  oscillating 
shear  flow;  Heat  transfer  case,  T  /T  =  1 


00*0 


Figure  8.  Comparison  of  preliminary  computations  (solid  line)  with  experimental 
results  for  mean  and  fluctuating  temperature,  and  for  the  temperature 
velocity  correlation  coefficients.  Molecular  Pr  =  0.7,  zero  pressure 
gradient . 
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